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ABSTRACT 


Increasing the transfer (HT) coefficient used in thermal industries is very important. Various methods are used 
to improve the efficiency of thermal heat HT so that maximum HT takes place in a smaller space. Ethylene 
glycol (EG) is generally used as an agent for convective HT. EG obtains energy from ahot source and discharges 
it to the required location. At present, the most consumption of EG is to produce engine cooling fluid. In the 
upcoming research, the TB of EG fluid in two-dimensional microchannels(MCs) has been investigated using 
molecular dynamics (MD) simulations, and the effect of variables such as MC dimensions and MC wall 
temperature(Temp) on the TB of the simulated fluid has been investigated. The results revealed that by 
increasing the Temp dif- ference of the MC wall from 10 to 50 K, the maximum temperature (Max-Temp) 
and velocity (Max-Vel) of the target sample increased to 640.94 K and 0.024 A/ps. It can be concluded that the 
increase in the cross-sectional area and the wall Temp difference leads to an increase in the HT rate in the 


MC. 


1. Introduction 


Nowadays, the problem ofimproving HT in engineering sciences and 
various industrial applications has been increasing. It has attracted the 
attention of many researchers, so it has become an essential part of 
experimental and theoretical research [1-3]. Increasing the HT coeffi- 
cientusedin thermalindustries is veryimportant. Various methods are 
used toincrease thermal HT efficiency so that maximum HT takes place 
ina smaller space. Improving HT using conventional methods has 
significantly saved energy costs and resources and preserved the envi- 
ronment [4-6]. The convection method mainly uses EG as an agent for 
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HT. EG is usually used in air conditioning, cold water systems, or cars 
with liquid cooling fluid. At present, the most consumption of EG is to 
produce engine cooling fluid. EG is used as an antifreeze in cars and 
machines. In the cooling system of cars, EG liquid plays the basic 
function of dealing with freezing. But due to the corrosiveness of this 
substance, anti-corrosion inhibitors are added to it to prevent rusting 
[7=10]: 

Many papers have been done on TB and improving the thermal 
performance of various fluids, and some have been briefly discussed. For 
instance, Guo et al. |11] investigated external heat’s effect on the base 
fluid’s TB in an MC. In this study, EG was considered as the base fluid. 


The obtained results revealed that an increase in Temp leads to an in- 
crease in entropy and subsequently the improvement of HT in the 
studied fluid. Bagheri et al. [12] investigated the forced convection HT 
in a thermal MC with H20 /alumina and H20 /copper oxide NFs. They 
discussed the effect of hydraulic diameter and volume fraction of NPs on 
HT. This study revealed that increasing the volume fraction of NPs and 
decreasing the diameter of NPs increases the value of the Nusselt 
number. Wu et al. [13] examined the effect of nanochannel roughness 
on the thermal performance of fluid in a nanochannel. The results 
revealed that increasing the nanochannel roughness can disturb the TB 


of the fluid. Kalteh et al. [14] investigated the effect of Nps’ volume 


Table 1 
The parameters of 1. LJ potential functions [27,28]. 


Particle type o (A) € (kcal/mol) 
H 2.886 0.044 

oO 3.500 0.06 

Cc 3.851 0.105 

Pt 2.754 0.080 


percentage and diameter on the TB of the studied NF in an MC. The 
obtained results revealed that the TB of NF improves by increasing the 
volume fraction of NPs and decreasing the diameter of NPs. Rostami 
et al. [15] investigated the effect of cubic roughness on the TB of argon 
fluid inside a platinum MC. Hu etal. [16] examined the effect of MC type 
on the atomic AB of fluid in an MC. The results revealed that changing 


the channel type can increase the atomic stability of the whole system. 
Shangetal.{17|examined the MCthicknessontheABofNF. Theresults 


revealed that with the increase of wall thickness, the atomic displace- 


ment decreases and the results indicate the aggregation process in the 


structure. Chen and Ding [18] analyzed an MC heat sink’s HT charac- 


teristics and cooling performance with H20 /alumina NFs with different 
vol magidiacd Qi TNE Sued that the effect of channel wall Temp, di- 
mensions, and cross-sectional area on the atomic and thermal perfor- 
mance of EG in a two-dimensional MC has not been investigated. 
Therefore, the TB ofthe desired fluid in the two-dimensional MC has 
been investigated in the current study. By studying and reviewing the 
research done sofar, itcan be seenthatfewstudies have beendonein 


determining the thermal properties of EG fluid in a two-dimensional 
platinum MC using the MD simulation {16-19}. Therefore, in this 
research, theeffects of MC wallTempdifferenceandMCdimensions and 


cross-sectional area on physical parameters suchas density, Temp, ve- 
locity, and thermal parameters suchas HF and TC have been evaluated 
using LAMMPS software and MD simulation. 


2. Simulation method 
2.1. The MD method 


Nowadays, MD is widely used in various engineering and basic sci- 
ences. Analytical determination of complex molecular systems’ char- 
acteristics is impossible due to a large number of particles [19,20]. 
Therefore, it can be solved using the calculation method defined in the 
MD simulation. Newton’s equations of motion are fundamental re- 
lationships in the MD. Solving these equations for each particle of the 
system leads todetermining the movement path ofall particles over time 
[21,22]: 

dU 
F,;=ma;= —V;U = — (1) 
dr; 

In relation (1), mjrepresents the mass of thei particle, and in the 
relations discussed, the acceleration is also obtained according to rela- 
tion (2) [21,22]: 
gate (2) 

' de 

In relation (2), r;is a vector that specifies the position of the particle. 
From the combination ofthe above relations, relation (3) is obtained 
[21,22]: 


— (3) 
dU dri 
Td de 


One ofthe standard numerical integration methods is the velocity- 


relations (4) and (5){23;24]: 
oa D aa 

r (t+ An = r (t4 At + (An)* +... (4) 
C C 
ra- A)= r A Mrt 2 TR (Ans... (5) 

The above relations, after simplification, can be expressed according 
to the following relation [23,24% > 

&r; 2 

r (t+ AD =2r(4)— r(t— Ad 4 "(Ad (6) 
i i i dt2 
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TC coefficient (K) via Green-Kubo relation [25,26]: 
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Inthe above equation, the parameterseiandVrespectivelyrepresentthe 
total energy (for each particle) and the total volume of the sample. Also, 
in this equation, Si represents the entropy of the system and vi represents 
the velocity of the particle. Another factor that is very important in the 
MD method is the definition of the potential function to determine the 
force between the constituent particles ofa structure [27]. 

Eq. (10) represents the potential function of Lennard Jones (LJ), 
which is used in the MD simulations to simulate interactions between 
different particles. This function has the potential of repulsion and 
attraction in its formulation [28]. 


() [¢ 12 C del 
Uri =4e Z — 2 sre (10) 


The parameters ofthe LJ potential function are revealed in Table 1. 
Also, based on formulas (11) and (12), the values of g and € are 


obtained [28,29]: 


Ej = Ys (11) 
o- GG (12) 
ij 2 


In EAM, force-fiejd atomic interaction is achieved as follows [30]: 


=) 12 ©) 
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U= Fa Pp Vij Pop Vij (13) 
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Verlet’s algorithm [23,24]. This integration algorithm is according to 
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Fig. 1. Schematic of a two-dimensional platinum MC from (a) side and (b) front view. 


250 


150 


Temperature {K} 


100 


‘Lime (ns) 


Fig. 2. Temp variations vs. the simulation time 


1.1. Simulation details 


This research used Avogadro and LAMMPS software to model the 
structure of EG and two-dimensional platinum MC. The dimensions of 
the SBin this study are 0.5.0.5 5d umŝ. The total number of present 
particles is 10,032 atoms. Itis worth noting that the Z-degree boundary 
conditions are considered fixed, while these conditions are considered 
periodicin both the X and Y directions. In the first stages, itis necessary 
tomake sure about the atomic modeling and the force field. In this step, 
the NVT ensemble (aconstant volume and Temp ensemble) is used. The 
simulation time and time step used in these studies equals 20 ns and 1 fs. 
For this purpose, the desired fluidis balanced atthe desiredinitialTemp 
(300 K). The initial Temp in the structures is balanced using the nose- 
Hoover thermostat. From the computational point of view, the 
balancing process in the current research was carried out for 10 ns. Fig. 1 
demonstrates a schematic of the modeled structure. In the next step, the 
effect of MC walldimensions and Temp on EG fluid’s TBin the simulated 
MCsisexamined. Inthis step, the NVEensemble(aconstantvolumeand 
energy ensemble) is used. Thesesteps were examined during 10ns. 


reported. Fig. 2 reveals the Temp changes in the defined atomic sample 
to the simulation time. The results for this quantity indicate a decrease in 
1.2. Stability process 


The changes in thermodynamic quantities are investigated and 


Temp fluctuations in this structure over time, indicating the 
appropriate settings applied in the current research. The Temp 
value in the whole atomic sample is equal to 300 K in the last time 
step. The Temp balance in the samples indicates the non-divergence 
of the oscillation amplitude of the structures within the SB, which 
indicates the stability of the structure. As a result, it can be said 
that the simulation time of 10 ns is sufficient in the balancing 
stage. 

Also, the KE changes reveal similar behavior, and this 
quantity converges to the value of 17.35 eVafter 10 ns. The changes 
inKEinthe atomicsample with simulated time are presentedin Fig. 3. 
Considering that the value of the second power of the velocity in the 
equation of KE is positive, as a result, this value is always positive. 
Convergence in KE with time results from decreasing mobility in 
atomic samples. By reducing the mobility of atoms inside the 
number of oscillations created decreases, and as a result, the structure 
reveals physical balance [31,32]. 


1.3. Validation 


In the final part of this part, it will be useful to provide a correct 
quantity to validate the samples. For this purpose, the yge of the simu; 
lated fluid radial distribution function (RDF) is presented in Fig.4. The 
use of RDF indicates how the atoms ina structure are arranged relative 
to each other. In the analysis of the result obtained in this part, itcan be 
said that there is a distinct peak and several times the liquid phase is 
revealedin the simulated base fluid sample, which somehowreveals the 
accuracy of the simulation method. 


2. Results and discussion 


As mentionedin the previous sections, the atomicand TB of EG fluid 
inside a two-dimensional platinum MC are investigated in the current 
research. For this purpose, in this section, the influence of factors such as 
the dimensions ofthe MC and the Temp difference of the walls on the 
atomicand TBofthe fluidis reported. Also, inthe present study, various 
physical quantities such as density, TC, velocity, and Temp and HF 
profiles have been calculated and reported. 


2.1. Effect of wall Temp difference 


After observing the simulated samples’ thermodynamic 
equilibrium, the sample’s atomic and TB are investigated. Creatinga 


Temp difference in the MC walls can effectively move fluid atoms near 
the MC walls and, 


finally, their atomic and TB. To investigate this issue, the Temp 
differ- ence between the two walls is considered equal to 10 K, 20 K, 
30 K, and 50 K. The thermophysical parameters of density, Temp, 
velocity, HF, and TC have been calculated and reported. Fig. 5 
reveals the atomic arrangement of EG fluid after applying a Temp 
difference of 10, 20, 30, and 50 K. According to Fig. 5, the 
simulated fluid is still stable by applying the Temp difference in 
the MC walls. 


Fig. 6 reveals the maximum density (Max-Dens) changes in the 
simulated samples according to the Temp difference between the MC 
walls. According to Fig. 6, as the Temp difference in the 
sampleincreases, the maximum value of the density increases, and this 


increase occurs near the cold wall. As the wall Te difference, 


increases and due to the higher Temp near the hot wall, the density 
decreases. The exis- tence of HT from the hot wall to the fluid 
causes the KE of the fluid 


molecules to increase and the time of their presence in the force field of 
the wall to decrease significantly. As a result, the density in the vicinity 
of the hot wall decreases. As the Temp difference between the walls 
increases from 0 to 50 K, the Max-Dens value decreases from 0.025 
atom/A3 to 0.3018 atom/ A. 


Fig. 7 illustrates the Max-Temp changes in the simulated samples 
according to the Temp difference between the MC walls. By increasing 
the wall Temp difference, the effect of the wall Temp on the distribution 
of atoms near the wall increases. As a result, the process of HTs between 
the walls increases compared to the fluid atoms in the vicinity of the 
wall. As a result, the HTs between the wall and the fluid increase. As a 
result, by increasing the Temp difference between the walls from 0 to 50 
K, the Max-Temp value increases from 514.35 K to 640.94 K. 

Fig. 8. represents the Max-Vel changes in the simulated samples ac- 
cording to the Temp difference between the MC walls. According to the 
presented Fig. 8, with the increase of the Temp difference between the 
wall from 0 to 50 K, the Max-Vel value increases from 0.01 A/psto 0.024 
A/ps. This increase is the result of increasing the HF applied to the 
particles, and as a result, their mobility and velocity increase. On the 
samples due to the increase in the Temp difference and the range of 
fluctuations. As a result, the Temp of the defined samples increases. This 
AB of samples should be considered in designing atomic structures for 
HT in industrial and confirmed applications. 

Afterobservingthe atomic balance andstabilityand determining the 
structural evolution in the samples, changes in HF and TC are investi- 
gated according to the Temp difference. Fig. 9 reveals the HF of the 
simulated samples. Based on these results, it can be said that the increase 
in Temp difference causes an increase in the oscillation range and, as a 
result, an increase in the HF flowing inside the SB. The value of this 
quantity at the Temp difference of 10, 20, 30, and 50 Kis equal to 1356, 
1388, 1403, and 1459 W/m? after 10 ns. 
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Fig. 3. KE variations vs. the simulation time. 
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Fig. 4. The simulated fluid RDF. 


Fig. 5. The atomic arrangement of EG fluid after applying Temp difference with values of (a) 10, (b) 20, (c) 30, and (d) 50 K. 
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Fig. 6. Max-Dens variations in the simulated samples with the Temp difference. 


651 
£ 
= 60d 
= 554 
= soa 
x 
2 


450 


Temperature Difference (K) 


Fig. 7. Max-Temp variations with the Temp difference. 
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Fig. 8. Max-Vel variations with the Temp difference. 


Fig. 10 in this section illustrates the TC of the simulated 
samples versus time in various Temp differences. By maximizing 
the HF by 


other hand, thestated process causes anincreaseinthemovement ofthe 
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Fig. 9. Changes in the fluid’s HF according to the MC wall’s Temp difference. 
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Fig. 10. Changes in the fluid’s TC based on the MC wall’s Temp difference. 


Table 2 


Numerical results of the simulated structure for the Temp difference of cold and 
hot walls. 


Temp Max-Dens Max-Vel Max-Temp HF (W/ TC (W/ 
difference (K)  (atom/A’) (A/ps) (K) m°?) m.K) 
10 0.023 0.013 561.359 1356 0.849 
20 0.021 0.018 592.011 1388 0.853 
30 0.020 0.021 621.065 1403 0.886 
50 0.018 0.024 640.943 1459 0.903 


applying the Temp difference of the MC wall, the TC also increases so 
thatthe value ofthe TC at the Temp differenceof10,20,30,and 50K 
after 10ns, respectively, is equalto0.849,0.853,0.886and0.903 W/m. 
K. As aresult, the HT in industrial cases is expected to increase by 
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increasing the Temp difference in the designed sample. 

The numerical results obtained in this part of the research are pre- 
sented in Table 2. Table 2 demonstrates the values of density, velocity, 
Temp, and thermal parameters (HF and TC) with increasing Temp dif- 
ferences of MC walls. 


2.2. The effect of MC dimensions 


The size of the cross-section of the defined MC is very influential in 
the fluid flow inside it. As aresult, itcan affect the simulated assembly’s 
atomic transformation and thermal conduction. To investigate this issue, 
the platinum MC’s cross-section size equals 0.25, 0.36, 0.49, and 0.56 
A?. Fig. 1 lrepresentstheatomicarrangementofEG fluidin different MC 
cross-sections. According to the atomic structures presented in Fig. 11, 
the fluid sample has stability with the increase of the cross-sectional area 


Fig. 11. The atomic arrangement of EG fluid according to the cross-sectional area of the MC (a) 0.25 A?, (b) 0.36 A?, (c) 0.49 A?, (d) 0.56 A?. 
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Fig. 12. Max-Dens changes according to MC cross-section. 


ofthe MC. Max-Dens decreases with the increase of the cross-sectional area. This 
Fig. 12 demonstrates the Max-Dens changes in the simulated samples behavior indicates that as the cross-sectional area of the MC increases, 

according to the cross-sectional area of the MC. Based on the results the flow intensity of EG decreases compared to the cross-sectional area 

obtained for the Max-Dens in the samples, it can be concluded that the attheTempof300K.Byincreasingthecross-sectionalareafrom0.2 5A? 
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Fig. 13. Max-Temp changes according to MC cross-section. 
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Fig. 14. Max-Vel changes according to MC cross-section. 


to 0.56 A”, the Max-Dens value reaches from 0.025 atom/ A3to0.016 
atom/A3. 

Increasing the cross-sectional area of the MC causes the absorbent 
force to enter the fluid particles with less intensity from the side of the 
walls. This increases the mobility of fluid particles, and as a result, the 
Max-Temp in the samples increases. Fig. 13 illustrates the Max-Temp 
changes in the simulated samples according to the cross-sectional area 
of the MC. According to the graph drawn in Fig. 13, with the increase of 
the cross-sectional area from 0.25 Å? to 0.56 A”, the Max-Temp value 
increases from 514.353 K to 691.639 K. This maximum value occurred 
in the middle regions of the MC, which had a greater distance from the 
walls. Because less forceis applied to the fluid particlesin this part from 
the side of the walls, they have the maximum movement and velocity. 

Therefore, the Temp has also reached the maximum value in this part. 

Similar behavior can be seen for the Max-Velin Fig. 14. Fig. 14 
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represents the Max-Vel changes in the simulated samples according to 
thecross-sectionalarea ofthe MC. Byincreasing the cross-section from 
0.25 A? to 0.56 Å?, the Max-Vel value increases from 0.01 A/psto0.014 

A/ps. The process described in this part demonstrates that by increasing 
thecross-sectional area, the power ofHTinthesamplecan beincreased, 
which can be expressed and proven by performing calculations related 
to the HF and TC of the simulated base fluid. 

In this part’s final stage, the samples’ HF and TC are calculated and 
reported according to the cross-sectional area of the platinum MC. 
Fig. 15 represents the changes in the HF of the EG fluid inside the 
platinum MC according to the cross-section of the MC. According to the 
results obtained from the graph, the sufficient time for establishing 
thermal equilibrium in the simulated sample is 10 ns. based on these 
graphs, itcan be said that 10 ns is enough time to establish thermal 
equilibrium in the simulated sample. Therefore after this time, the 
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Fig. 15. Changes in the HF according to the cross-sectional area of the MC. 
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Fig. 16. Changes in the TC according to the cross-sectional area of the MC. 


Table 3 


Numerical results of simulated structure in terms of MC cross-section. 


MC cross Max-Dens Max-Vel Max-Temp HF (W/ TC (W/ 
section (A?) (atom/A*) (A/ps) (K) m°’) m.K) 
0.25 0.025 0.010 514.353 1253 0.845 
0.36 0.022 0.012 589.331 1286 0.848 
0.49 0.018 0.013 618.001 1349 0.859 
0.56 0.016 0.014 691.639 1381 0.868 


amount of HF in the sample with a cross-sectional area of 0.25, 0.36, 
0.49, and 0.56 Å? are equal to 1253, 1286, 1349, and 1381 W/m?, 
respectively. Asaresult, with theincreaseofthecross-sectionoftheMC, 
the amount of HF inside it increases. This increase is defined as the result 
ofanincreasein theatomic mobility and the oscillating amplitude ofthe 
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fluid. 

Considering that there is a direct relationship between the trans- 
ferred HF and the conductivity coefficient ofatomic samples, the in- 
crease ofthe HF willcorrespondtotheTC. Fig. 1 illustrates the changes 
in the TC of the EG fluid inside the platinum MC according to the cross- 
section of the MC. The value of the TC in the sample has an increasing 
trend with the increase ofthe cross-section of the channel. The numer- 
ical results represent an increase in the TC from 0.845 to 0.868 W/m.K 
W/m.K due to the increase in the cross-sectional area of the MC from 
0.25 to 0.56 A? (see Table 3). 

The Numerical results of the platinum MC and EG fluid sample in 
terms ofMC cross-sectional areaat the Temp of 300 Kandafter 10nsare 
presented in Table 3. 

Generally, the MCs can be used in the field of pumps, microchip 
cooling, air conditioning, etc. Therefore, it is possible to improve the TB 


inside these MCs by adjusting the wall Temp or channel dimensions. 
3. Conclusion 


The current research investigated EG fluid’s atomic and TB inside the 
platinum MC. For this purpose, computer simulations were used, and the 
MD simulation method was used to be more precise. From a technical 
point of view, the simulations carried out in this research included two 
main stages balancing the atomic samples and establishing the atomic 
and thermal transformation. The results obtained in the balancing stage 
indicated the thermodynamic balance of the structure in the defined 
initial conditions (Temp of 300 K). The equilibrium process was inves- 
tigated in 10 ns. In the second step, the samples were studied again for 
10 ns, and their atomic and TB was checked using the Green Kobo 
method. The overall results obtained in this part ofthe research were as 
follows: 


. Balance in atomic samples due to the convergence of physical 
quantities, including Temp and KE, is 300K, 17.35 eV. On the other 
hand, these convergences indicate that the duration of 10 ns is suf- 
ficient to balance the structures. 

. Increasing the Temp difference between the wall from 0 to 50 K has 
increased the Max-Temp and Max-Vel from 514.35 to 640.94 K 
(about 25%) and from 0.01 to 0.024 A/ps (more than twice), 
respectively. 

. Byincreasing the Temp difference between the wall from 10 to 50K, 
the amount of HF and TCincreased from 1356 W/m?and0.849 W/m. 
K to 1459 W/m? (about 8%) and 0.903 W/m.K (about 7%). As a 
result, the HT in industrial cases is expected to increase by increasing 
the Temp difference in the designed sample 

. By increasing the cross-sectional area from 0.25 A” to 0.56 A’, the 
Max-Temp and Max-Vel value increased from 514.353 to 691.639 
(about 34%) K and from 0.01 to 0.014 A/ps (about 40%). So, by 


incr asing the crass-sectional area, the power of HT in the sample 
can be irtereased. 


- Increasing the cross-sectional area from 0.25 A to 0.56 & has 
increased the values of HF and TC from 1253 W/m énd 0.845 W/m. 


K to 1381 W/m? (about 10%) and 0.868 W/m.K (about 3%). This 
increase is defined as the result ofan increase in the atomic mobility 
and the oscillating amplitude of the fluid. 
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